Chemistry in Id via DMRG: Benchmarks for small atoms, ions, and molecules 
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The density matrix renormalization group (DMRG) method, an efficient solver of lattice models 
in Id, has been adapted to solve real-space problems. We use DMRG to study the behavior of 
soft-Coulomb interacting matter, reporting accurate benchmarks for small systems. We compare 
with Hartree-Fock and the local density approximation, noting similarities and differences with 3d. 
We study how correlation grows as bonds stretch. 
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I. INTRODUCTION AND PHILOSOPHY 

Electronic structure methods such as density func- 
tional theory (DFT) are excellent tools for investigat- 
ing the properties of solids and molecules — except when 
they are not. Standard density functional approxima- 
tions in the Kohn-Sham (KS) framework [1] work well 
in the weakly correlated regime [2-4] , but these same ap- 
proximations can fail miserably when the electrons be- 
come strongly correlated [5]. Though DFT is exact in 
principle, what can these successes and failures tell us 
about the fundamental nature of KS-DFT? In particu- 
lar, do these failures indicate that approximate KS-DFT 
cannot describe strongly correlated systems? 

Strongly correlated electrons are usually studied with 
lattice Hamiltonians, which directly model electron local- 
ization and spin interactions. Both are correlation effects 
arising from an underlying continuum model but which 
are not readily apparent — nor easy to treat — within most 
methods applicable to real space. As simplified mod- 
els of bulk materials, lattice Hamiltonians offer qualita- 
tive insight into generic behavior, such as phase tran- 
sitions and collective states of matter. They are often 
sufficient for extracting universal features, such as scal- 
ing exponents near phase transitions, but do not yield 
highly accurate estimates of material-specific properties, 
such as bond lengths and atomization energies. Among 
many clever methods for solving lattice models, DMRG 
(density-matrix renormalization group) is outstanding 
for its accuracy and speed, especially for one-dimensional 
problems, and has been applied to highly-entangled sys- 
tems with hundreds of lattice sites [6, 7]. 

An intriguing approach to the problem of strong corre- 
lation has been the development of density functional ap- 
proximations for lattice models. A DFT for lattice mod- 
els, sometimes called site-occupation functional theory 
(SOFT), can be formulated and applied to these prob- 
lems [8]. The most basic energy functional, the local 
density approximation (LDA) , can be found from certain 
limits solvable by the Bethe ansatz [9]. SOFT approxi- 
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mations for more exotic quantities such as the entangle- 
ment entropy have been tested [10, 11], and exact time- 
dependent Kohn-Sham DFT has been investigated [12]. 
Progress in this area provides new approximate methods 
for solving lattice models at a fraction of the cost of tra- 
ditional many-body approaches, but with all the usual 
uncertainties associated with approximate DFT calcula- 
tions. The main use of such calculations is really as a 
test of functional ideas on systems (albeit idealized) that 
have strong correlation. 

However, simplified lattice models are not good mimics 
of realistic systems for the properties that are extracted 
from typical DFT calculations. Take for example a two- 
site Hubbard model of the hydrogen molecule. Using a 
minimal basis of a Is orbital on each H atom, this model 
becomes relatively exact as the separation tends to infin- 
ity, and recovers the Heitler-London wavefunction. Near 
equilibrium separations, this model becomes empirical, 
as the parameters are not uniquely determined. 

The two-site Hubbard model also highlights a criti- 
cal difference between the nature of DFT approximations 
and SOFT approximations. For this model, the SOFT 
energy functional becomes a function of two numbers (the 
site occupations) whose sum is constrained by the total 
particle number. Thus it becomes a simple function. Yet 
most of the subtleties (and excitements) of DFT approx- 
imations are because the Hilbert space is infinite, and we 
use functionals, not functions. Theorems can be proven 
on lattices that are not true for the exact continuum func- 
tionals [13]. Many of the challenges of modern functional 
approximations concern the derealization error [14] (a 
special case of which is self-interaction error [15]). The 
nature of this error is completely masked when entire 
atoms are described by single sites. Furthermore, we can 
rationalize the success of density functional approxima- 
tions by considering a certain, physically-important scal- 
ing of the coordinates: in the semiclassical limit of all 
continuum systems, LDA becomes exact [16, 17]. This 
limit can never be attained by a site model with finite- 
range interactions. It appears that, for SOFT approxi- 
mations to succeed, they must do so for fundamentally 
different reasons than real-space DFT approximations. 

The work described here is a crucial step toward an 
entirely new approach to the pressing problem of dealing 
with strong correlation, with an eye toward bulk materi- 
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FIG. 1. The KS potential for a stretched hydrogen molecule 
found from interacting electrons in Id. 



als. Instead of attempting to learn about functionals by 
defining density functions on lattice models, we adapt an 
extremely powerful lattice-model solver to work for the 
continuum. We have built a library of routines that al- 
low application of DMRG to Id systems with real-space 
potentials and long-range interactions, instead of lattice 
models. Compared to past applications of DMRG, this 
setting far more closely mimics three-dimensional real- 
ity, and especially those aspects that are captured by 
present popular DFT approximations, i.e., the semiclas- 
sical limit. In the present work, we present accurate 
benchmark calculations for small systems, both approxi- 
mate and exact. These establish the qualitative similar- 
ities and differences between this Id world and the 3d 
world, and demonstrate how basic questions concerning 
the performance of density functionals can be directly 
addressed with this method. 

Our results are illustrated in Fig. 1, which shows Id H2 
with soft-Coulomb interactions. The exact density was 
found by DMRG and inverted to find the correspond- 
ing exact KS potential. The bond has been stretched 
beyond the Coulson-Fischer point, and a strong XC con- 
tribution to the KS potential is needed to reproduce the 
exact density in the bond region [18]. Such calculations 
have often been performed for few-electron systems in 3d 
in the past [19, 20], but our method allows exact treat- 
ment of much larger systems. In another paper [21], we 
show how powerful this method is, by solving a chain 
of 100 Id H atoms, and by performing interacting inver- 
sions. All such calculations were previously unthinkable 
for systems of this size, and unreachable by any other 
method. We have applied these techniques to perform 
the first ever Kohn-Sham calculations using the exact XC 
functional, essentially implementing the exact Levy-Lieb 
constrained search definition of the functional, which we 
will present in yet another paper. 



II. BACKGROUND IN DMRG 

The density matrix renormalization group (DMRG) is 
a powerful numerical method for computing essentially 
exact many-body ground-state wavefunctions [6], Tra- 
ditionally, DMRG has been applied to Id and quasi-2d 
finite-range lattice models for strongly correlated elec- 
trons [7]. DMRG has also been applied to systems in 
quantum chemistry, where the long-range Coulomb inter- 
action is distinctive. The Hamiltonians which have been 
studied in this context include the Pariser-Parr-Pople 
model [22] and the second-quantized form of the Hartree- 
Fock equations, where lattice sites represent electronic 
orbitals [23, 24]. 

DMRG works by truncating the exponentially large 
basis of the full Hilbert space down to a much smaller- 
one which is nevertheless able to represent the ground- 
state wavefunction accurately. Such a truncation would 
be highly inefficient in a real-space, momentum-space, or 
orbital basis; rather, the most efficient basis consists of 
the eigenstates of the reduced density matrix computed 
across bipartitions of the system [6]. A DMRG calcu- 
lation proceeds back and forth through a Id system in 
a sweeping pattern, first optimizing the ground-state in 
the current basis then computing an improved basis for 
the next step. By increasing the number of basis states 
m that are kept, DMRG can find the wavefunction to 
arbitrary accuracy. 

The computational cost of DMRG scales as N s m 3 
where N s is the number of lattice sites. For gapped sys- 
tems in Id, the number of states m required to compute 
the ground-state to a specified accuracy is independent 
of system size, allowing DMRG to scale linearly with N s . 
For gapless or critical systems, the m needed grows log- 
arithmically with system size, making the scaling only 
slightly worse. The systems considered here have a rela- 
tively low total number of electrons such that the number 
of states m required is small, often less than 100. This in 
turn enables us to work with the very large numbers of 
sites (N s ~ 1000 — 5000) needed to reach the continuum, 
as described in more detail below. 



III. METHODOLOGY: CHEMISTRY IN ID 

To apply DFT in its natural context — in the 
continuum — we shall consider a model of soft-Coulomb 
interacting matter [25-27], where the electron repulsion 
has the form 

«cc(w) = 7===p (1) 

Vtt + 1 

and the interaction between an electron and a nucleus 
with charge Z and location X is 

v{x) = -Zv ee {x-X). (2) 

The soft-Coulomb interaction is chosen to avoid diver- 
gences when particles are close to one another, and 
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has been used to study molecules in intense laser fields 
[25, 26]. The wavefunctions and densities within this 
model lack the cusps present in 3d Coulomb systems. 
However, the challenge presented by the long-range in- 
teractions in 3d Coulomb systems remains for these Id 
model systems. 

Most DMRG implementations for solving Id systems 
have been limited to lattice models with finite-range in- 
teractions. Yet we would like to study continuum systems 
with long-range interactions, systems such as those en- 
countered in quantum chemistry. Although DMRG has 
been applied in quantum chemistry by having lattice sites 
represent basis functions [23, 24], such an approach re- 
quires a basis appropriate to a particular molecule and is 
not flexible enough to represent an arbitrary real-space 
system. 

Instead, we enable DMRG to operate in the continuum 
by discretizing over a fine real-space grid. With a lattice 
spacing of a, the real-space Hamiltonian for a Id system 
becomes in second quantized notation, 

+ vJ n i + \ Y v ™ n ' 1 ( n J ~ ' ( 3 ) 

3 ij 

where fx = ii — 1/a 2 , v J = v(j a) and v l J c = v cc (\i — j\ a). 
The S{j in the last term cancels self interactions. The 
operator cj CT creates (and Cj a annihilates) an electron of 
spin a on site j, fij = + riji, and rij a — cj^c^. 

The hopping terms cj^Cj+i^ (and complex conjugate) 
come about from a finite-difference approximation to the 
second derivative. Like the second-quantized Hamilto- 
nians considered in quantum chemistry, this Hamilto- 
nian corresponds to an extended Hubbard model; Eq. 
(3), however, is motivated from a desire to study the Id 
continuum alongside familiar DFT approximations. Be- 
cause we require that the potentials and interactions vary 
slowly on the scale of the grid spacing, the low-energy 
eigenstates of the discrete Hamiltonian (3) will approxi- 
mate the continuum system to very high accuracy. More- 
over, we check convergence with respect to lattice spac- 
ing. Because our potentials — and thus our ground-state 
densities — vary slowly on the scale of the grid spacing, we 
can accelerate convergence by using a higher-order finite- 
difference approximation to the kinetic energy operator; 
this simply amounts to including more hopping terms in 
Eq. (3)._ 

Even in its discretized form the Hamiltonian Eq. (3) 
represents a challenge for DMRG because of the long- 
range interactions. Including all 7V S 2 interaction terms, 
where N s is the number of lattice sites, would make the 
calculation time scale as iV s 3 overall. Fortunately, an el- 
egant solution has been recently developed [28] which 
involves rewriting the Hamiltonian as a matrix prod- 
uct operator (MPO) a string of operator-valued matri- 
ces. This form of the Hamiltonian is very convenient 
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c — » oo 


0.1000 


-81.50 


-82.30 


-82.40 


-82.41 


0.0500 


-19.58 


-20.46 


-20.57 


-20.58 


0.0200 


-2.22 


-3.16 


-3.27 


-3.29 


0.0100 


0.27 


-0.68 


-0.80 


-0.82 


0.0050 


0.90 


-0.07 


-0.18 


-0.20 


0.0025 


1.06 


0.09 


-0.03 


-0.05 


-> 


1.12 


0.14 


0.02 


0.00 



TABLE I. Convergence of model hydrogen energy with re- 
spect to lattice spacing a and distance c from the atom to 
the edge of the system, with differences in units of micro- 
hartree from the infinite continuum extrapolation of E = 
-0.66977714. 

for DMRG, and MPOs naturally encode exponentially- 
decaying long-range interactions [29] . Assuming that our 
interaction v ee (u) can be approximated by a sum of expo- 
nentials, the the calculation time scales only linearly with 
the number of exponents -ZV exp used. This reduces the 
computational cost from 7V S 3 to N s N cxp . In practice, for 
our soft-Coulomb interactions and modest system sizes 
(N s < 1000), we find that only -/Vc X p = 20 exponentials 
are needed to obtain an accuracy of 10~ 5 in our approxi- 
mate v ee (u). The largest -/V exp we use in this paper is 60, 
which is necessary to find the equilibrium bond length of 
Id H2 accurate to ±0.01 bohr (a system with N s « 2000). 

For technical reasons, we take all of our systems to have 
open (or box) boundary conditions. This has no adverse 
effect on our results because we can extend the grid well 
past our edge atoms. The extra grid sites cost almost no 
extra simulation time due to the very low density of elec- 
trons in the edge regions. To evaluate the dependence of 
the energy on these edge effects and the grid size, con- 
sider Table I. This table shows the convergence of the 
Id model hydrogen atom ground-state energy with re- 
spect to the lattice spacing a and the distance c from the 
atom to the edge of the system, using the second-order 
finite difference approximation for the kinetic energy, as 
in Eq. (3). Our best estimate for the Id H atom energy is 
-0.66977714, converged to at least microhartree accuracy, 
which differs slightly from that of Eberly [25] . 

In addition to the accurate many-body solutions of- 
fered by DMRG, we can also look at approximate so- 
lutions given by standard quantum chemistry tools. 
Hartree-Fock (HF) theory can be formulated for these 
Id systems by trivially changing out the Coulomb in- 
teraction for the soft-Coulomb. The exchange energy is 
then: 

E x = - i Y Y J dx J dx' v cc (x - x') x 

a i,j = l 

0ia(x)<Pja (x)(f)ja (x'). (4) 

In performing HF calculations, instead of using an orbital 
basis of Gaussians or some other set of functions, our 
"basis set" will be the grid, as in Eq. (3). This simple and 
brute force approach allows us a great degree of flexibility, 
but is only computationally tractable in Id. 
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In this setting we also implement DFT. As mentioned 
in the introduction, DFT has been applied directly to 
lattice models. But our model and interaction are meant 
to mimic the usual application of DFT to the contin- 
uum. In particular, the LDA functionals we will use are 
similar to their 3d counterparts, unlike the Bethe ansatz 
LDA (BALDA) , which has a gap built in [9] . One calcu- 
lates the LDA exchange energy by taking the exchange 
energy density per electron for a uniform gas of density 
n, namely e" mf ( n )j an d then integrating it along with the 
electronic density: 



} = J dxn(x)el nil (n(x)). 



(5) 



We find e™ if (n) by evaluating Eq. (4) with the KS or- 
bitals of a uniform gas. For a uniform gas, the KS or- 
bitals are the eigenfunctions of a particle in a box, whose 
edges are pushed to infinity while the bulk density is 
kept fixed. Because the interaction has a length-scale, 
i.e. v cc (ju) 7 p u co (u) for some p, even exchange is not 
a simple function. One finds: 



^unif 



(n) = -n/(fc P )/2, 



where k F — 7m/ 2 is the Fermi wavevector and 

1 



/(*) = 



dy 
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sin y 



V 2 v/z 2 + y 2 ' 
In fact, / is related to the Meijer G function 



/(*) = G 
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/(4s) 



(6) 



(7) 



(8) 



We write r s = l/(2n) as the average spacing between 
electrons in Id. In Fig. 2, we show the exchange energy 
per electron for the unpolarized g RS elS cl function of r s . 
For small r s (high density), e™ if -> -1/2 + 0.203 r s ; for 
large r s (low density), e™ if -> -0.291/r s - ln(r s )/(4r s ). 
For contrast, in 3d, the exchange energy per electron is 
always -0.458/r s [31], where r s = (3/(47rn)) 1 / 3 . 

In practice, we do not use pure DFT, but rather spin- 
DFT, in which all quantities are considered functionals 
of the up and down spin densities. In that case, we 
need LSD, the local spin-density approximation. For ex- 
change, there is a simple spin scaling relation that tells 
us [32] 



,unif 



(n t ,n+) = -n (1 + a() 2 /(fc F (l + aC))/4, (9) 

CT=±1 

where £ = (n-f — n^)/n is the polarization. This is less 
trivial than for simple Coulomb repulsion. At high den- 
sities, there is no increase in exchange energy due to spin 
polarization, while there is a huge increase (tending to a 
factor of 2) at low density, as shown by the solid black 
line in Fig. 2. In fact, e™ if (r s , C = 1) = e£ nif (r s /2, C = 0). 

To complete LDA, we need the correlation energy den- 
sity of the uniform gas at various densities and polariza- 
tions. We are very fortunate to be able to make use of 



the pioneering work of Ref. [33], which performs just such 
a QMC calculation and parametrizes the results, yield- 
ing accurate values for e" nlf (r s ,£), which are also plot- 
ted for the unpolarized and fully polarized cases in Fig. 
2. These curves are not qualitatively similar to the 3d 
e"c lf ( r s:0- For these Id model systems, the fully polar- 
ized electrons almost completely avoid one another at the 
exchange level, so that correlation barely decreases their 
energy for any value of r s . For unpolarized electrons, the 
effect of correlation is to make them avoid each other en- 
tirely for low densities (r s > 5) and the XC energy per 
electron becomes independent of polarization. However, 
for unpolarized electrons at high density, correlation van- 
ishes with r s , and exchange dominates, as in the usual 
3d case. For moderate r s values, the correlation contri- 
bution grows with r s , as shown by the red dashed line of 
Fig. 2. To give an idea of what range of r s is important, 
for the hydrogen atom of Fig. 3, 95% of the density has 
r s (x) — (2n(x))~ 1 between 1 and 8. 
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FIG. 2. Parametrization of the LDA exchange and exchange- 
correlation energy densities per electron for polarized f = 1 
and unpolarized ( — densities [33]. 

Armed with these parametrizations and tools, we are 
ready to discover Id quantum chemistry. 



IV. RESULTS 

DMRG gives us an excellent tool for finding exact an- 
swers within a model Id chemistry world. Our Id world is 
designed to mimic qualitatively the 3d world, not match 
it exactly. Below we explain some important differences 
between our model Id systems and real 3d systems, start- 
ing with the simplest element. 



A. One-electron atoms and ions 

We begin with Id hydrogen and calculate its ground- 
state energy to high precision. We find that E(H) — 
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FIG. 3. The hydrogen atom with both exact and LSD den- 
sities, as well as the LSD KS potential. 

—0.66977714, accurate to 1 microhartree. Table I shows 
the convergence with respect to the lattice spacing and 
the size of the system. Eberly et al. [25] were the first to 
consider the soft-Coulomb hydrogen atom and its eigen- 
states. Its ground-state energy is similar to the 3d hy- 
drogen atom energy of —0.5 a.u. Because the potential 
and wavefunction is much smoother, the kinetic energy is 
only 0.11 a.u., as opposed to 0.5 a.u. in 3d. Because the 
potential does not scale homogeneously, in Id the virial 
theorem does not yield a simple relation among energy 
components, unlike 3d. 

Again because of the lack of simple scaling, hydrogenic 
energies do not scale quadratically for our system. A 
simple fit of energies for Z > 1 yields: 

E z « -Z + VZ/2-2/9 + a 1 /Vz, (N = 1) (10) 

where otx = 0.0524 is chosen to make the result accurate 
for Z = 1. The first two coefficients are exact in the large- 
Z limit, where the wavefunction is a Gaussian centered 
on the nucleus. 

A well-known deficiency of approximate density func- 
tionals is their self-interaction error. Because E x is ap- 
proximated, usually in some local or semilocal form, it 
fails to cancel the Hartree energy for all one-electron sys- 
tems. Thus, within LSD, the electron incorrectly repels 
itself. This error can be quantified by looking at how 
close £ , ^ SD is to the true E x . As can be seen in Table 
II, E^ SD is about 10% too small. For hydrogen, the self- 
interaction error is about 30 millihartrees. By adding in 
correlation, this error is slightly reduced, but remains fi- 
nite. This is an example of the typical cancellation of 
errors between exchange and correlation in LSD. 

As a result of self-interaction error, the LSD electron 
density spreads out too much, as shown in Fig. 3. In this 
figure we can also see how the LSD KS potential fails 
to replicate the true KS potential, which for hydrogen is 



system 


T 


E E h 


Ex E x E c 


H 
He+ 
Li++ 
Be+++ 


0.111 
0.192 
0.258 
0.316 


-0.670 -0.647 
-1.483 -1.455 
-2.336 -2.304 
-3.209 -3.176 


-0.346 -0.311 -0.007 
-0.380 -0.343 -0.006 
-0.397 -0.359 -0.005 
-0.408 -0.369 -0.005 



TABLE II. Exact and LSD results for Id one-electron sys- 
tems. 



the same as the external v(x). And although the LSD 
KS potential is almost parallel to v(x) where there is a 
large amount of density, it decays too rapidly as \x\ — > oo. 
What this adds up to, both in Id and in 3d, is that LSD 
will not bind another electron easily, if at all. We will 
return to this point when considering anions. 

B. Two-electron atoms and ions 

For two or more electrons, the HF approximation is 
not exact. The traditional quantum chemistry definition 
of correlation is the error made by HF: 

E$ c = E - E HF . (11) 

In Table III, we give accurate energy components for two- 
electron systems; recall that the components do not sat- 
isfy a virial theorem in our Id systems. The total energy 
can be fit just as for one-electron systems, but now: 

E Z tt-2Z + yfZ + CQ-a 2 /VZ, (N = 2) (12) 

where cq = 0.507 and a-i = 0.235. The HF energies 
may be fit with c^ F = 0.476 and af F = 0.167. These 
fits are not accurate enough to give the large Z behav- 
ior of E^ c , which seems to vanish as Z — > oo. For 3d 
two-electron systems, the correlation energy scales to a 
constant at large Z [34]. Overall, |£^? C | is much smaller 
in Id than in 3d. Rather than the dimensionality, it is 
the soft nature of our Coulomb interactions that causes 
the reduction in correlation energy compared to 3d. The 
exact wavefunctions in 3d have cusps whenever two elec- 
trons of opposite spin come together, caused by the di- 
vergence of the electron-electron interaction. This cusp- 
related correlation is sometimes called dynamic correla- 
tion; any other correlation, involving larger separations 
of electrons, is called static [35]. (Note that the dis- 
tinction between static and dynamic correlation is not 
precise.) Our soft-Coulomb potential has no divergence 
and induces no cusps, so dynamical correlation is min- 
imal. There is little static correlation in tightly bound 
closed shell systems, such as our Id Li + and Be ++ , so 
\Ec | <C \E\. In contrast, for H~, where one electron 
is loosely bound, one expects most of the correlation to 
be static even in 3d, and one sees large and similar E c 
values in Id and 3d. In Section IV E, we discuss some 
quantitative measures of strong correlation. 

Next we study the exact Kohn-Sham DFT energy com- 
ponents of these two-electron systems. Here we need the 
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system 


T V Vee 


E E Hb E£ 


H" 
He 
Li+ 
Be ++ 


0.115 -1.326 0.481 
0.290 -3.219 0.691 
0.433 -5.084 0.755 
0.556 -6.961 0.790 


-0.731 -0.692 -0.039 
-2.238 -2.224 -0.014 
-3.896 -3.888 -0.008 
-5.615 -5.609 -0.006 


3d H~ 
3d He 
3d Li+ 
3d Be++ 


0.528 -1.367 0.311 
2.904 -6.753 0.946 
7.280 -16.13 1.573 
13.66 -29.50 2.191 


-0.528 -0.488 -0.042 
-2.904 -2.862 -0.042 
-7.280 -7.236 -0.043 
-13.66 -13.61 -0.044 



TABLE III. Exact and HF two-electron atoms and ions, in 
1- and 3-d (exact data from Ref. [19], Li + is fit quadratically 
to surrounding elements, and HF data from Ref. [36, 37]). 





Ts U £ xc 


E x E c T c 


H" 
He 
Li+ 
Be++ 


0.087 1.103 -0.595 
0.277 1.436 -0.733 
0.425 1.542 -0.779 
0.551 1.601 -0.806 


-0.552 -0.043 0.028 
-0.718 -0.014 0.013 
-0.771 -0.008 0.008 
-0.801 -0.006 0.005 


3d H" 
3d He 
3d Li+ 
3d Be++ 


0.500 0.762 -0.423 
2.867 2.049 -1.067 
7.238 3.313 -1.699 
13.61 4.553 -2.321 


-0.381 -0.042 0.028 
-1.025 -0.042 0.037 
-1.656 -0.043 0.041 
-2.277 -0.044 0.041 



TABLE IV. Energies of the exact KS system for two-electron 
atoms and ions. 3d data (Li + fitted) from Ref. [19]. 



DFT definition of correlation, which differs slightly from 
the traditional quantum chemistry version: 

E c = E - (T s + V + U + E x ) 

= T C + U C , (13) 

where E x is the exchange energy of the exact KS or- 
bitals, T s is their kinetic energy, U is the Hartree en- 
ergy, T c = T — T s is the kinetic correlation energy, and 
Uc = V cc — U — E x is the potential correlation energy. 
All these functionals are evaluated on the exact ground- 
state density, with numerical results found in Table IV. 
The difference between the quantum chemistry Ec and 
the DFT E c is never negative and typically much smaller 
than \E C \ [38]. For the two-electron systems of Table III 
and Table IV, the difference is zero to the given accuracy 
for all atoms and ions besides Id H - . For our systems, 
just as in 3d, Ec C — E c vanishes as Z — > oo. All the large 
DFT components (T s , U, E xc ) are typically smaller than 
their 3d counterparts and scale much more weakly with 
Z. However, our numerical results suggest T c —E c as 
Z — > oo, just as in 3d. 

To obtain the KS energies for a given problem, we 
require the KS potential, which is found by inverting 
the KS equation. For one- or two-electron systems, this 
yields: 

v s (x) = —^*^{x). (N<2) (14) 
2 v / n(x) ax 1 

For illustration, consider the exact KS potential of Id 
helium in Fig. 4. Inverting a density to find the KS po- 
tential has also been done for small systems in 3d, where 







Id He / \ 


i 

— exact 
LDA 








v!r DA (x) 

- v s [n](x) 
- v(x) 
i 





-4 -2 2 4 



x 

FIG. 4. The exact KS potential for a model helium density 
found from interacting electrons in Id, as well as the LDA 
density and LDA KS potential found self-consistently. 



QMC results for a correlated electron density have proven 
extremely useful [19]. One can find simple and useful 
constraints on the KS potential by studying the large 
and small r behavior of the exact result [19, 39]. In 3d, 
for large r, the Hartree potential screens the nuclear po- 
tential, and the exchange-correlation potential goes like 
— 1/r [40]. In Id, the softness of the Coulomb potential 
is irrelevant, so the Hartree potential screens the nuclear 
potential for large |x| as in 3d. Though it seems likely, 
we have no proof that the exchange-correlation potential 
for the soft-Coulomb interaction should tend to — l/|x| 
for large \x\, analogous to the 3d Coulomb result. To 
check this would require extreme numerical precision in 
the density far from the atom, due to the need to eval- 
uate Eq. (14) where the density is exponentially small. 
Instead, in Figures 1 and 4, we require the KS poten- 
tial to go as — 6/ 1 a; | once the density becomes too small 
(around n.« 10~ 5 , which happens at |x| ~ 6 for helium), 
and we fit b based on what the potential (14) is doing 
before that point. The actual value of b has no visible 
effect on the density on the scale of these figures. 

We now consider the performance of LDA for these 
two-electron systems, starting with how well LDA repli- 
cates the true KS potential. Though the LDA density is 
only slightly different from the exact density on the scale 
of Fig. 4, the LDA potential clearly decays too rapidly 
(exponentially) at large r and is too shallow overall, just 
as in 3d [19]. Like the hydrogen atom discussed ear- 
lier, this is a result of self-interaction error. LDA energy 
results are given in Table V. Clearly LDA becomes rela- 
tively more accurate as Z grows, because XC becomes an 
ever smaller fraction of the total energy. Comparing Ta- 
bles IV and V, we also see that LDA underestimates the 
true X contribution by about 5%, while overestimating 
the correlation contribution, so that XC itself has lower 
error than either, i.e., a cancellation of errors. 
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E LDA % error 


^xc ^x i^c 


ET 
He 
Li+ 
Be ++ 


-0.708 -3.1% 
-2.201 -1.7% 
-3.850 -1.2% 
-5.564 -0.9% 


-0.601 -0.536 -0.065 
-0.690 -0.646 -0.044 
-0.731 -0.696 -0.035 
-0.753 -0.723 -0.030 


3d H" 
3d He 
3d Li+ 
3d Be++ 


-0.511 -3.2% 
-2.835 -2.4% 
-7.143 -1.9% 
-13.44 -1.2% 


-0.419 -0.345 -0.074 
-0.973 -0.862 -0.111 
-1.531 -1.396 -0.134 
-2.082 -1.931 -0.150 



TABLE V. LDA for 2 electron systems. H~ does not con- 
verge in Id or 3d; the results are taken from using the LDA 
functional on the HF density [41]. 3d LDA data from Engel's 
OEP code [42]. 



x 

T3 



T3 




FIG. 5. The differential logarithmic decay of the helium atom 
density for various methods. The horizontal, dashed lines 
correspond to the asymptotic decay constants. 



Much insight into density functionals has been gained 
by studying the asymptotic decay of densities and poten- 
tials far from the nucleus [39]. In Fig. 5, we plot dn/dx/n 
to emphasize the asymptotic decay of the exact, LDA, 
and HF helium densities. The HF density is very accu- 
rate compared to the LDA density. For large x, the HF 
density has very nearly the same behavior as the exact 
density, and both approach their asymptotes very slowly. 
By contrast, the LDA density reaches its asymptote by 
x 4. For each approximate calculation, its asymptotic 
decay constant 7 can be found using the highest occupied 
molecular orbital (HOMO) energy: 7 = — 2enoMO- 
The asymptote for the exact curve can be found using 
the ionization potential / — E(N — 1) — E(N) of the sys- 
tem, which determines the density decay: 7 = —2y/2I. 
Because the HF asymptote lies nearly on top of the ex- 
act asymptote, Koopmans' theorem — or / = — c^omo — 
is extremely accurate for Id helium. We list both HOMO 
and total energy differences in Table VI. 

There is a long history of studying two-electron ions 
in DFT, including the smallest anion H _ , which presents 
interesting conundrums for approximate functionals [43, 
44] . Looking at the ionization energies of these 2 electron 





LSD 


HF 


exact 


system 


-ehomo I 


-ehomo I 


/ =-£HOMO 


H" 
He 
Li+ 
Be++ 


0.062 
0.478 0.747 
1.242 1.546 
2.064 2.389 


0.054 0.022 
0.750 0.741 
1.557 1.552 
2.404 2.400 


0.061 
0.755 
1.560 
2.406 



TABLE VI. Id HOMO eigenvalues and ionization potentials 
of two-electron atoms and ions, for the exact functional, LDA, 
and HF. LDA does not converge for H~ anion, but LDA en- 
ergies can be found using the HF density [41]. 



systems, we can extrapolate the critical nuclear charge 
necessary for binding two electrons, i.e., figuring out the 
Z value for which 1 = 0. This happens around Z = 0.90 
in Id, and around Z = 0.91 in 3d [45]. Within LDA, the 
critical value is above Z = 1, because H _ will not bind. 
DFT approximations have a hard time binding anions — 
both in Id and in 3d — due to self-interaction error. A way 
to circumvent this problem is to take the HF anion, which 
binds an extra electron, and evaluate the LDA functional 
on its density. As seen in Table VI, this approach is far 
better than either taking total energy differences or the 
negative of the HOMO energy from HF alone, just as in 



3d [41]. As in 3d, — Chomo ^ s useless as an approximation 
1 1 and J HF are close to each 



c HOM 

to /. The HF results -e 



HOMO 



other and closest to / for larger Z; but / LSD does very 
well for small Z, and is best for H _ . 



C. Many-electron atoms 

Before looking at larger atoms, a word of caution. In 
3d systems, degeneracies in orbitals with different an- 
gular momentum quantum numbers produce interesting 
shell structure. In Id, there is no angular momentum, 
and thus a lack of the shell structure we are familiar 
with. This means that it is not entirely clear which real 
elements our model Id atoms correspond to. The first 
three Id elements might well be called hydrogen, helium, 
and lithium; but the fourth Id element may behave more 
like neon than beryllium; to be consistent with Ref . [33] , 
we call it beryllium. To showcase Id Be, consider its 
LDA treatment in Fig. 6. The exact KS potential is also 
plotted, and the LDA KS potential roughly differs only 
by a constant in the high density region, just as with hy- 
drogen (Fig. 3) and helium (Fig. 4). In the low density 
regions, the LDA correlation potential is the dominant 
piece of the LDA KS potential. 

As we increase the number N of electrons in our sys- 
tems, correlation also increases, but HF theory is still 
better than LDA until N = 4. Exact and HF data for 
many-electron atoms can be found in Table VII, and LDA 
data in Table VIII. Despite good agreement with all other 
data, we did not find He - or Li~ to bind as in Ref. [33], 
neither in HF nor DMRG, nor LSD (not surprisingly). 
When using energy differences to calculate the ionization 
energy, HF outperforms LSD until beryllium, as can be 
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- f(x) 
(j) 2 (x) 

■— v s (x)/2 
v(x)/4 
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x 

FIG. 6. An LDA "beryllium" atom, complete with the LDA 
orbitals 4>j(x), LDA KS potential v^ T>A (x), exact KS poten- 
tial Vs(x), and external potential v(x). The density n LDA (a;) 
was found self-consistently using the LDA method, and barely 
differs from the true n(x) on this scale. 



system 


T V Ke 


E E Hb ' E^ 


Li 

Be+ 
Be 


0.625 -6.484 1.648 
0.922 -9.240 1.864 
1.127 -11.13 3.219 


-4.211 -4.196 -0.015 
-6.454 -6.445 -0.010 
-6.785 -6.740 -0.046 



TABLE VII. Exact and HF many-electron atoms and ions, 
in Id. 



seen in Table IX. For these systems, 7 LSD > I > J HF : 
LSD overestimates the ionization energy, and HF under- 
estimates it — just as in 3d. As with the fewer electron 
systems, the LSD HOMO energies are not a good way to 
estimate /, whereas the HF HOMO energies are. 





E LSD % error 


77LSD irnLSD r LSD 
-C-XC E x E c 


Li 

Be+ 
Be 


-4.179 -0.8% 
-6.410 -0.7% 
-6.764 -0.3% 


-1.044 -1.004 -0.041 
-1.117 -1.086 -0.031 
-1.450 -1.376 -0.075 



TABLE VIII. LSD energies for many-electron Id systems. 

To find the KS energy components for these many- 
electron (N > 2) systems, we again require the exact 
KS potential. For these systems, Eq. (14) is no longer 
valid, so we must find the KS potential another way. The 
simplest procedure is to use guess-and-check, adjusting 
the KS potential until its density can no longer be dis- 
tinguished from the target density found using DMRG. 
Updates to the KS potential can be more or less sophis- 
ticated, without changing the final result in the region 
where the density is large; in the low-density region, how- 
ever, two very different KS potentials can give rise to 
densities that are indistinguishable on the scale of our 
figures. However, the KS energy components do not rely 





LSD 


HF 


exact 


system 


-ehomo I 


-ehomo I 


I =-£homo 


Li 


0.166 0.329 


0.316 0.308 


0.315 


Be+ 


0.628 0.846 


0.842 0.835 


0.839 


Be 


0.162 0.353 


0.313 0.295 


0.331 



TABLE IX. Many-electron ionization energies for LSD, HF, 
and exact Id systems. 





T s U £ xc 


E x E c T c 


Li 
Be+ 
Be 


0.611 2.749 -1.087 
0.912 3.042 -1.168 
1.091 4.736 -1.481 


-1.071 -0.016 0.014 
-1.157 -0.011 0.009 
-1.430 -0.051 0.036 



TABLE X. Energies of the exact KS system for many-electron 
Id atoms and ions. 



significantly on the behavior of the KS potential out in 
the low-density region. For example, noticeable irregu- 
larities occur for the beryllium KS potential, but these 
occur beyond \x\ = 7 and are not plotted in Fig. 6. In 
Table X, the exact KS energies for some many-electron 
systems are tabulated. For Li and Be + , spin-DFT is 
used, but the spin-dependent energy components (such 
as TIf ) are summed together to give a spin- independent 
energy. 

The study of the energies of neutral atoms as N = 
Z — >• oo is important due to the semiclassical result being 
exact in that limit [16]. In this limit, the oldest of all 
density functional approximations, Thomas-Fermi (TF) 
theory, becomes exact [46]. However, due to a lack of 
scaling within the soft-Coulomb interaction, the large Z 
limit of the energy is non-trivial, making a semiclassical 
treatment difficult. A plot of the neutral atom energies 
as a function of N appears in Fig. 7. On this scale, both 
the LDA and HF results lie nearly on top of the exact 
curve. 



D. Equilibrium properties of small molecules 

We now briefly discuss small molecules near their equi- 
librium separation. In order to find the equilibrium bond 
length for our Id systems, we take the nuclei to be in- 
teracting via the soft-Coulomb interaction, just like the 
electrons. Given this interaction, consider the simplest of 
all molecules: the H^ cation. HF yields the exact answer, 
and LSD suffers from self-interaction (more generally, a 
delocalization error [5]). A plot of the binding energy is 
found in Fig. 8. Because the nuclear-nuclear repulsion 
is softened, the binding energy does not diverge as the 
internuclear separation R goes to zero. As seen in Table 
XI, LSD overbinds slightly and produces bonds that are 
too long between H atoms, which is also the case in 3d. 
The curvature of the LSD binding energy is too small 
near equilibrium, which makes for inaccurate vibrational 
energies, especially in 3d. This can also be seen in Table 
XL Finally, we note that the energy of stretched H^ does 
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N 

FIG. 7. Energies of neutral atoms in Id. 
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FIG. 8. The binding energy curve for our Id model ILJ , shown 
with an absolute energy scale, and with nuclear separation R; 
horizontal dashed lines indicate the energy of a single H atom. 



not tend to that of H within LSD, due to derealization 
error [5]. 

Next we consider H2. A plot of the binding energy is 
found in Fig. 9; the large R behavior will be discussed in 
the following section. Just as in 3d, HF underbinds while 
LDA overbinds; HF bonds are too short, and LDA bonds 
are too long. Further, HF yields vibrational frequencies 
which are too high, and LDA are a little small, which is 
the case both in Id and 3d. All of these properties can 
be seen in Table XL 



E. Quantifying correlation 

It is often said that DFT works well for weakly cor- 
related systems, but fails when correlation is too strong. 
Strong static correlation, which occurs when molecules 
are pulled apart, is also identified with strong correlation 
in solids [5]. Functionals that can accurately deal with 





HF 


LSD 


exact 


system 


D c (eV) 


T-r-L- 

HJ 


3.88 (0%) 


4.00 (3%) 


3.88 


3d 


2.77 (0%) 


2.89 (4%) 


2.77 


H 2 


2.36 (-23%) 


3.53 (15%) 


3.07 


3d H 2 


3.54 (-25%) 


4.80 (1%) 


4.75 


system 


Ro 




2.18 (0%) 


2.28 (4%) 


2.18 


3d 


2.00 (0%) 


2.18 (9%) 


2.00 


H 2 


1.50 (-6%) 


1.63 (2%) 


1.60 


3d H 2 


1.41 (1%) 


1.47 (5%) 


1.40 


system 


uj (xlO 3 cm" 1 ) 




2.2 (0%) 


2.0 (-9%) 


2.2 


3d Hj 


2.4 (0%) 


1.9 (-21%) 


2.4 


H 2 


3.3 (6%) 


3.0 (-3%) 


3.1 


3d H 2 


4.6 (5%) 


4.2 (-5%) 


4.4 



TABLE XL Electronic well depth D c , equilibrium bond ra- 
dius -Ro, an d vibrational frequency u) for the ILj~ and H 2 
molecules, with percentage error in parentheses. Exact 3d 
H 2 results taken from Ref. [47]; the remaining 3d values are 
from Ref. [36] using the aug-cc-pVDZ basis set [48]. 
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FIG. 9. The binding energy curve for our Id model H 2 , 
shown on an absolute energy scale, with nuclear separation 
R. Dashed curves represent unrestricted calculations. 

strong static correlation in stretched molecules can also 
accurately yield the band gap of a solid [49, 50]. Most 
DFT methods, however, fail in these situations. To see 
these effects in Id, we shall now examine three descrip- 
tors of strong correlation, which will be when no cor- 
relation is present and close to 1 when strong correlation 
is present. 

A simple descriptor of strong correlation is simply to 
calculate the ratio of correlation to exchange: 



In the limit of weak electron-electron repulsion, a goes to 
zero for closed-shell systems, and HF becomes exact. For 
example, for the two-electron atoms and ions in Table IV, 
a goes to zero as Z increases. In Table XII, we compute a 
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for various bond lengths of the hydrogen molecule, both 
in Id and in 3d. At the equilibrium bond length, a is 
small, indicating that the HF solution is very close to the 
exact. When the bond is stretched to R — 5, a increases 
ten-fold: a standard HF solution for a bond length of 
R = 5 does not do well at all. The Id and 3d results are 
remarkably similar. We can also compute a using the 
LDA functionals for E c and E K evaluated on the LDA 
density; however, a LDA is not as good of an indicator for 
strong correlation as the true a is. 

The second descriptor of strong correlation requires 
first understanding where correlation comes from. From 
Eq. (13), correlation can be separated into two pieces: (1) 
the kinetic correlation energy T c = T — T s , due to the 
small difference between the true kinetic energy and the 
KS kinetic energy, and (2) potential correlation energy, 
U c = V cc — U — E x . In the limit of weak correlation in 
3d, U c -» -2T C , so the ratio: 



Eg+Tg 

E a 



(16) 



has always been found to be positive, and vanishes in the 
weakly correlated limit [51], which we have also observed 
in Id. But if T c -C \E C \, we have correlation without 
the usual kinetic contribution, which occurs when sys- 
tems have strong static correlation. For example, in the 
infinitely stretched limit of H2 , T s — > T while E c remains 
finite, so /3 — > 1. In Table XII, we see that j3 increases as 
we stretch the H 2 molecule, both in Id and 3d. Thus f3 
is a natural measure of static correlation in chemistry. 

There is another test for strong static correlation, 
which we can use on closed-shell systems. In the ab- 
sence of a magnetic field, the ground-state wavefunction 
will be a spin singlet. Yet when stretching a molecule 
such as H 2 , there is a bond length beyond which an ap- 
proximate calculation, such as HF and LDA, will pre- 
fer to break spin symmetry in order to lower the en- 
ergy. This separation is known as the Coulson-Fischer 
point [52]. Beyond this point, an unrestricted calcula- 
tion, whose up orbitals may be different than down or- 
bitals, spontaneously breaks spin symmetry. In this case, 
we say the system has strong static correlation, because 
the electrons want to localize on different atoms. Be- 
fore the Coulson-Fischer point, however, both restricted 
and unrestricted calculations agree, and the system does 
not have strong static correlation. This phenomenon is 
observed in Fig. 9, where the solid (dashed) curves repre- 



sent restricted (unrestricted) calculations for a stretched 
hydrogen molecule. Unrestricted HF/LSD breaks spin 
symmetry around R = 2.1/3.4, beyond which the unre- 
stricted solution gives accurate total energies, but very 
incorrect spin densities [53]. The numbers are similar in 



1 = 2.3/3.3) 


[541. 












Id 


3d 




R 


1.6 


3.4 


5.0 


1.4 5.0 


exact 


a 


0.04 


0.21 


0.46 


0.06 0.45 




P 


0.21 


0.58 


0.87 


0.18 0.89 


LDA 


a 


0.09 


0.16 


0.21 





LABLE XII. Lable of correlation descriptors a and /3 
(Eqs. (15) and (16)) for H2 at an equilibrium and a stretched 
bond length R. 3d data from Ref. [55]. 



V. CONCLUSION 

In this paper we have set benchmarks for our Id model 
chemistry set. Though this Id model world is not the 
same as our 3d world, many properties are similar. Our 
DMRG results demonstrate that standard approxima- 
tions like HF and LDA behave just as in the real world. 
But our approach allows us to study strong correlation 
that arises directly from the continuum Hamiltonian. 
Understanding electron correlation in this Id world will 
give us insight into real 3d systems, and will illuminate 
the challenges that we face in making approximate DFT 
work for strongly correlated systems. Other approaches, 
such as range-separated hybrids [56], dynamical mean- 
field theory [57], and LDA+U [58] will be tested in the 
future. 

Using these tools, we seek to develop partition density 
functional theory (PDFT) [59]. PDFT is a promising 
method to break up molecules into atoms in a quan- 
titative and exact way. Much like the KS potential, 
which makes a set of non-interacting electrons look like 
they are interacting, a new partition potential turns non- 
interacting fragments (such as atoms) into interacting 
fragments. Hitherto, PDFT has mostly been studied 
with non-interacting electrons. With the machinery we 
have developed using DMRG, this will change. Through 
studying the exact theory, we hope to understand exact 
and useful constraints on the partition potential and the 
corresponding partition energy. 

With gratefulness, we acknowledge DOE grant DE- 
FG02-08ER46496 (KB and LW) and NSF grant DMR- 
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